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An approximate solution to the time-dependent density functional theory (TDDFT) response 
equations for finite systems is developed, yielding corrections to the single-pole approximation. 
These explain why allowed Kohn-Sham transition frequencies and oscillator strengths are usually 
good approximations to the true values, and why sometimes they are not. The approximation 
yields simple expressions for Gorling-Levy perturbation theory results, and a method for estimating 
expectation values of the unknown exchange-correlation kernel. 
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Traditional density-functional theory (DFT) is a pop- 
ular and efficient method for the calculation of ground- 
state properties of interacting many electron systems. 
While it has long been the method of choice for solid- 
state calculations recent improvements in approxi- 
mate functionals have also made it popular in quantum 
chemistry, because of the ability to handle large systems 
1^. Within exact ground-state theory, the Kohn-Sham 
eigenvalues and orbitals are mathematical constructs, 
designed to reproduce the ground-state density. The 
one exception is that the highest occupied eigenvalue is 
known to equal minus the ionization potential of the sys- 
tem. For many years, the unoccupied Kohn-Sham levels 
of solids have been used to interpret the excitations in 
solids, despite the well-known underestimate of the band 
gap of insulators in LDA. 

There exist several extensions of the basic formalism 
that allow extraction of excited-state properties. One 
popular approach is via time- dependent DFT (TDDFT), 
in which the interacting system in a time-dependent ex- 
ternal field is mapped exactly to a non-interacting time- 
dependent Kohn-Sham (KS) system with the same time- 
dependent density j|] . If a weak time-dependent electri- 
cal field is considered, this leads to a Dyson-like response 
equation for the exact susceptibility of the interacting 
electronic system. Poles of this susceptibility occur at 
the true transition frequencies and the strengths of 
these poles are related to oscillator strengths. Solution of 
these equations has been implemented in several quan- 
tum chemical packages, and excitation spectra of many 
molecules have been calculated and reported in hundreds 
of papers (See references in Ref. |^). Once an accurate 
ground-state Kohn-Sham potential is used, transition fre- 
quencies are typically within about 0.2 eV of experiment. 
Oscillator strengths are usually good in these calculations 
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(within a factor of 2), but not always. There are consider- 
able subtleties in applying TDDFT to extract the optical 
response of solids . 

The beauty of using TDDFT to extract optical spec- 
tra is that it combines moderate accuracy with inexpen- 
sive calculation, just as in the ground-state case. But 
even beginning from an exact ground-state KS potential, 
the spectrum still depends on the unknown exchange- 
correlation (XC) kernel, i.e., the functional derivative of 
the time-dependent XC potential. This is often approx- 
imated by the crude but reliable adiabatic local-density 
or generalized-gradient approximations. These TDDFT 
calculations compare very favorably with configuration- 
interaction singles (CIS) calculations, the only alterna- 
tive that is comparable in computational cost Q] . Higher 
level calculations, such as more complete CI, Bethe- 
Salpeter, or quantum Monte Carlo can be made more 
accurate, but cost more, limiting their use to smaller sys- 
tems. For example, the TDDFT approach has recently 
been applied to electron-transfer problems in biological 
systems 

Although TDDFT methodology has been implemented 
and is being used widely, understanding of its accuracy 
and reliability has been slow, as well as its relation to 
other methods. The relation to first-order Gorling-Levy 
perturbation theory has been shown [|o[ |l^, as well as 
the connection with the GW approximation The 
extreme case of stretched H2 has recently been studied by 
several authors 1 1^, |l^, Q , although this also represents 
difficulties for the ground-state theory jl^. By using 
a matrix formulation of Casida |l6| , the present paper 
shows how, when the excitations of a system are discrete, 
an approximate solution, that can be made arbitrarily 
accurate, can be used to understand and explain many 
trends in the results of TDDFT calculations. 

To demonstrate our results, we apply them to the pro- 
totype systems of the He and Be atoms. We chose these 
because their exact ground-state KS potentials are known 
p^t , and because their lowest allowed transitions exhibit 
two classes of systems we are interested in. In the He 
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atom the Is to 2p singlet transition is at 21.22 eV, while 
the KS transition, i.e., the energy difference between the 
Is and 2p ground-state KS orbital energy levels, is 21.15 
eV, less than 0.1 eV smaller. On the other hand, in the Be 
atom, the 2s to 2p level singlet transition is at 5.3 eV, but 
in the KS case it is at 3.6 eV. We explain below the fun- 
damental difference between these two systems, and why 
the KS eigenvalues are a good approximation in the first 
case, but a poor approximation in the second. Further- 
more, the oscillator strength for the Is to 2p transition in 
the He atom is 0.27, but 0.32 in the exact KS case. Thus 
the oscillator strength of the He atom KS system is close 
to the true one, but noticeably less close than the tran- 
sition frequency. Finally, the 2s to 2p oscillator strength 
in the Be atom is 1.37, but 2.54 in its Kohn-Sham alter 
ego, far closer than one would expect, given the error in 
transition frequencies. We explain these facts. 

We denote the exact KS eigenvalues as and orbitals 
as ^i(r). Casida |]l6| has written the TDDFT response 
equations as an eigenvalue equation for the square of the 
transition frequencies: 



(1) 



where g is a double index, representing a transition from 
occupied KS orbital i to unoccupied KS orbital a, ujq = 
— e^, = Lj^, and <i>g(r) = 0*(r)(/)a(r). The matrix is 



qq' 



Sqq' -|- 2 



t^q^'q{'l\fnxcii^)W) 



(2) 



where 



('?|/HXc(t^)k') = 



■'<f*(r)/Hxc(r,r',c.)$,,(r'). 

(3) 

In this equation, /hxc is the Hartree-exchange-correlation 
kernel, l/|r — r'| -|-/xc(r, r', w), where /xc is the unknown 
XC kernel. 

It has been noticed that the KS transition frequencies 
are often 'good' approximations to the true frequencies 
. If the transition frequencies and oscillator strengths 
are expanded in powers of the coupling constant A, as 
in GL perturbation theory ||l^ , the zero-order values are 
the Kohn-Sham values @. In Eq. (|), we see that, if 
/hxc is small, i.e., if the system is sufficiently weakly cor- 
related, this is correct. We show below that this is not 
the reason why the KS values are good approximations in 
many systems, such as the He atom, and is especially un- 
true when the ground-state has a near-degeneracy, such 
as for the Be atom. 

Our approximate solution relics on the fact that 
(q I /hxc I?') decays rapidly with distance from the diago- 
nal, because the overlap of increasingly different orbitals 
decay by cancellation of oscillations. To zero order, we 
ignore all off-diagonal elements, finding|^l| the small- 
matrix approximation (SMA)j2^: 



n^flq + 2LL>q{q\fiixc{^q)\q). 



(4) 



The original single-pole approximation (SPA) Q can be 
viewed as a special case of SMA when the shift from the 
KS value is small: 

= \J^q + 2Wg(g|/Hxc|'7) ^^q + (^l/nxclt^q)!?) + ■ • ■ 

(5) 

In the special case of including just the exchange ker- 
nel, this result is identical to first-order GL perturbation 
theory |l§. 

To go beyond SPA, we use a continued-fraction method 
p3| for inverting a matrix with a dominant diagonal. 
Truncating the expansion at second order in the off- 
diagonal matrix elements: 



SMA 



4c.,..,,|(g|/HXcK^^^)kOP 



QSMA _ fiSMA 



This is a key result of this paper, leading to many 
conclusions. First, it yields the exact GL perturbative 
expression to second-order in A. Expanding /hxc = 
A/hx + A2/1" + ..., wefind 

LU^LOq + X{q\f„^{u;q)\q) + X^Su;f\ (7) 

where the second-order shift consists of four terms: 

= (,i/g'K)i,)+2^ -^''^y-^;^^)'^'^'^ 



9' 7^9 

+ (g|/HxK)|g)(g|^^K)|g) - 



K'zI/hxK)!?)!^ 



2uJa 



(8) 



For the ground-state energy, the second-order correction 
has been identified as playing a key role in constructing 
accurate functionals, especially in cases of strong static 
correlation |2^. It is likely to play a similar role for 
excitations, and can be easily extracted from Eq. (|^). 

Second, we may now deduce precisely when SPA (or 
SMA) is valid. Defining the shift from the KS value as 



Aiiq = n-nq, 



(9) 



we rewrite Eq. (||) in the following suggestive form: 



SMA 



A^iSMA 



\M, 



qq' I 



QSMA _ ^SMA MqqMq.q, 



q'^q 1 



where M„, 



(10) 

(9 1 /hxc 1 9')- A simple estimate of 
the size of this correction can be given by assuming 
Mqq' ^ MqqMqiqi . This would bc cxactly true if 
/Hxc(r, r', 1.^) — f{r,uj)f{r',uj), and is accurate to within 
0.5% for exchange in the He atom. Then the SMA is 
valid when 
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ArjSMA 



flSMA _ QSMA 
q'^q 1 1 
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TABLE I: Exact results for the He and Be atoms, using 
numerically exact ground-state Kohn-Sham potentials. 



to leading order, we find: 



Atom 


trans. 


frequency 


(eV) 


osc. 


str. 






KS 


SMA" 


exact 


KS 


exact* 


He 


ls-^2p 


21.15 


21.23 


21.22 


0.3243 


0.2762 




Is— >3p 


23.06 


23.10 


23.09 


8.47(-2) 


7.34(-2) 




^4p 


23.73 


23.75 


23.75 


3.41(-2) 






Is— >5p 


24.04 


24.05 


24.05 


1.71(-2) 


1.50(-2) 




ls^6p 


24.21 


24.22 


24.22 


9.8(-3) 


8.6(-3) 


Be 


2s^2p 


3.61 


4.95 


5.28 


2.5422 


1.3750 




2s->3p 


7.33 


7.39 


7.46 


3.79(-2) 


9.01(-3) 




2s->4p 


8.29 


8.31 


8.33 


2.06(-2) 


2.3(-4) 




2s^5p 


8.69 


8.70 


8.69 


1.08(-2) 


8.1 (-4) 




2s^6p 


8.90 


8.90 


8.90 


6.3(-3) 


7.5(-4) 



"Hybrid SPA results from Ref. 
'He numbers are from Ref. lEf 



Q, converted to SMA. 
I and Be numbers from Ref. [ pvt . 



i.e., the SMA shift need only be small on the scale of 
the separation between transition frequencies. Thus, 
even when the corrections to KS transition frequencies 
are large, SMA can remain valid if the poles are well- 
separated. 

In the special case when the SMA correction is small 
compared to a KS transition itself, this result simplifies 
to 



{q\Uxc\q') 



< 1. 



(12) 



For example, the allowed transitions from the ground- 
state of the He atom are listed in Table |. Comparing KS 
transition frequencies with physical ones, we find them 
good to within less then 0.1 eV. This implies all matrix 
elements {qlfnxcW) are small, and SPA is valid. Calcu- 
lations within ALDA for this case show no difference 
between full solution of the response equations and the 
SPA result. The column marked SMA hsts SMA results 
with our best estimate of /hxc for this case, a hybrid of 
exact exchange with ALDA antiparallel correlation [pT| . 
For the transition to 2p, the exact exchange result is 21.37 
eV 1 28 , showing that there must be substantial cancel- 



lation by correlation effects. For this system and others 
like it, GL perturbation theory converges slowly, while 
our expansion converges rapidly. 

Our other prototype is the excitations from the ground 
state of the Be atom. For the 2s 2p transition, the ex- 
pectation value of /hxc is relatively large. We expect 
SMA to work quite well for that transition, but less well 
for others, since the one strong transition contributes to 
the correction in Eq. (10), especially having a small de- 
nominator. This is born out by the frequency results in 
Table |. Within ALDA the SMA transition is at 
5.27 eV, but the full calculation is 5.08 eV. 

In the SMA, in which off-diagonal elements are ne- 
glected, the eigenvalues in Eq. (|I]) remain unit vectors, 
and the oscillator strengths retain their KS values. When 
we include the change due to the off-diagonal elements 



E 



4(q|/Hxck')^g'^9'M«A*9 



f^SMA 
9 



^SMA 



(13) 



Here flq denotes the KS dipole matrix element. Our first 
conclusion from this important result is that it contains 
the exact GL expression for oscillator strengths to first 
order in A, by ignoring the correlation kernel. The sum is 
rapidly converging, as the matrix element decays rapidly 
with principal quantum number. Using /x for the He 
atom, the corrections of Eq. (p^), summed over only 
bound-bound transitions, reduce the oscillator strength 
by 11%, whereas the exact answer is 15% lower than 
the KS value. The remaining reduction is due to either 
correlation effects or transitions to the continuum. 

If we estimate off-diagonal elements with geometric 
means of diagonal elements, we find 



Al^SMA^^SMA 



1^ . (14) 



The effect of off-diagonal matrix elements is to mix vari- 
ous KS oscillator strengths. For the dominant transition 
, if Eq.(ll) is satisfied for excitation energies, it is also 
satisfied for oscillator strengths. The correction to an os- 
cillator strength is first-oidei in the off-diagonal matrix 
element, as opposed to the second-order correction to the 
SMA transition-frequency shift. Thus, fractional correc- 
tions to KS oscillator strengths will generally be larger 
than those to SMA shifts. So even with large SMA shifts, 
the associated oscillator strengths can be good. 

The oscillator strengths of the He and Be atoms con- 
firm our previous conclusions. For the well-separated 
transitions (He atom), the KS oscillator strengths are 
close to the true oscillator strengths, but not as close as 
the transition frequencies. The deviations estimated in 
Eq. ( jl^ ) are consistent with those of the transition fre- 
quencies. Similarly, in the case of the Be atom, we see 
that the 2s^2p KS oscillator strength is good to within 
a factor of 2, because the corrections due to other tran- 
sitions are quite mild. On the other hand, the higher 
transitions have KS oscillator strengths that are an or- 
der of magnitude different from the true ones, because of 
the huge corrections due to the first transition. 

Ours is not an accurate numerical solution of the 
TDDFT response equations, but is rather a method for 
understanding results. For example ||l^, the He oscil- 
lator strengths within ALDA are good to within 5%. 
The present work shows that this reflects the accuracy 
of (qlf^xc'^W) foi' these transitions. Our approximation 
scheme handles well-separated poles. As the system-size 
grows, e.g., even for small clusters, excitations become 
more dense on the w-axis, and kernel corrections can 
cause levels to reorder, etc. We believe that these meth- 
ods can still provide estimates of expected trends, but 
this remains to be tested. It is unclear what happens 
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in the thermodynamic hmit, without better approxima- 
tions to the XC kernel[^. For finite systems, the Tamm- 
DancoS^ approximation neglects all matrix elements be- 
tween one forward and one backward transition (reduc- 
ing the size of matrices by a factor of 2), and has been 
found to often yield excellent results compared to a full 
solution^. However, this approximation violates the 
Thomas-Reiche-Kuhn sum rule, whereas our expansion 
satisfies it order-by-order. 

We have focussed exclusively on the allowed transi- 
tions, because we need exact results on both transi- 
tion frequencies and oscillator strengths for our analy- 
sis. To analyze the limitations of SMA for singlet-triplet 
splittings would require knowledge of the exact oscil- 
lator strengths of spin-flipping transitions. 

So far, we have focussed on approximate solutions of 
the TDDFT response equations for the exact ground- 
state KS potential. In practice, this potential is approx- 
imated. Local-density and generalized-gradient approx- 
imations have potentials that are too shallow, so that 
Rydberg states are not bound. This can be corrected by 
some addition of the correct asymptotic behavio r [p9| , or 
by use of an orbital-dependent functional ||3^, |31|, p2[, 
whose derivative yields an accurate potential at large 
distances. From accurate calculations on the He and 



Be atoms, we find the principal effects of using either 
exchange-only or LDA-SIC potentials to be a shift in the 
orbital energies, numerically identical to the error in the 
ionization potential. But the KS dipole matrix elements 
are extremely accurate in these approximate potentials, 
so that the dominant error in oscillator strengths comes 
from the errors in eigenvalues. Thus current technology 
allows accurate calculation of KS oscillator strengths. 

We conclude with an observation that should be use- 
ful for development of approximations to /xc- In cases 
where poles are well-separated and SMA is valid, as can 
be determined by comparing oscillator strengths, the dif- 
ference between KS transition frequencies and physical 
ones yields an exact expectation value of /hxc, to all 
orders in coupling-constant A. This would provide an in- 
valuable benchmark for testing approximate XC kernels, 
similar to the widespread use exact XC potentials have 
enjoyed in the ground-state case[p7|. 
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